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Summary 

Two-scale  continuum  equations  are  derived  for  heterogeneous  continua  with  full  nonlinear 
electromechanical  coupling  using  nonlinear  mathematical  homogenization  theory.  The  resulting 
coarse-scale  electromechanical  continuum  equations  are  free  of  coarse-scale  constitutive 
equations.  The  unit  cell  (or  Representative  Volume  Element)  is  subjected  to  the  overall  mechanical 
and  electric  fields  extracted  from  the  solution  of  the  coarse-scale  problem  and  is  solved  for 
arbitrary  constitutive  equations  of  line-scale  constituents.  The  proposed  method  can  be  applied  to 
simulate  the  behavior  of  electroactive  materials  with  heterogeneous  line-scale  structure  and  design 
new  electroactive  materials  and  devices. 

Keywords:  electroactive,  multiscale,  homogenization,  electromechanical  coupling,  unit  cell, 
electostriction,  composite,  polycrystalline,  fine-scale  structure. 


1.  Introduction 

All  industrial  materials  are  heterogeneous  at  some  scale.  The  overall  behavior  of  materials  is 
controlled  by  their  fine-scale  structure  and  materials  with  tunable  line-scale  structure  can  exhibit 
superior  performance  at  coarse  level.  Multiphase  line-scale  structure  allows  improving  useful 
properties  while  reducing  the  effect  of  undesirable  properties  of  line-scale  constituents.  The  overall 
properties  can  be  tailored  by  changing  line-scale  structure  (mixing  different  materials  at  different 
scales,  changing  the  size  and  shape  of  constituents,  etc)  to  meet  the  needs  of  specific  applications. 
Attempts  have  even  been  made  to  state  an  optimization  problem  aimed  at  constructing  an  optimal 
material  architecture  [1,2, 3, 4, 5, 6, 7, 8].  Due  to  the  synergy  of  fine-scale  constituents  it  is  possible  to 
obtain  new  properties  which  are  not  available  in  fine-scale  constituents  alone.  Metamaterials 
[9,10,11,12,13]  represent  one  fascinating  example  of  how  tuning  fine-scale  structure  provides  a 
completely  different  dynamic  behavior  of  materials. 

Historically,  the  primary  goal  of  composite  materials  with  tunable  fine-scale  structure  was  to 
reduce  weight  in  particular  in  aerospace  applications.  More  recently,  emphasis  has  been  placed  on 
electroactive  materials,  which  found  their  use  in  actuators,  transducers  and  structures  capable  of 
changing  their  properties  (behavior)  depending  on  the  environment.  These  are  so-called  “smart 


structures”  where  feedback  loop  comprises  of  sensors  and  actuators.  The  presence  of  biasing  fields 
makes  such  materials  apparently  behave  differently  [14,15],  which  creates  an  opportunity  to 
control  the  response  of  so-called  smart  structures  to  various  excitations.  Electroactive  materials 
are  used  for  actuators,  sensors  and  smart  structures.  Applications  range  but  not  limited  to 
underwater  acoustics,  biomedical  imaging  [16],  robotic  manipulations,  artificial  muscles 
[17,18,19],  active  damping,  resonators  and  filters  for  frequency  control  and  selection  for 
telecommunication,  precise  timing  and  synchronization  14],  MEM  and  NEM  devices  [20,21,22], 
flow  control  [23],  precise  positioning  [24,  25],  conformal  control  surfaces  [26],  power  electronic 
devices  [27],  computer  memory  applications  [21,22],  tunable  optics,  generators  for  harvesting 
energy,  tactile  sensors  [28,29,30,31],  etc. 

Composites,  combining  electroactive  materials  with  other  materials,  allow  to  overcome  the 
limitations  intrinsic  to  pure  electroactive  materials  such  as  low  actuation  capabilities,  high  electric 
fields,  fast  decay  and  instabilities.  The  overall  properties  of  such  composites  are  extremely 
sensitive  to  fine-scale  details.  Various  studies  [19,26,32,33,34,35]  have  shown  that 
electromechanical  coupling  can  be  significantly  improved  by  making  non-homogeneous 
electromechanical  actuators,  composites  of  flexible  and  high  dielectric  modulus  or  even 
conductive  materials.  The  overall  response  of  a  composite  actuator  can  be  vastly  superior  to  that  of 
its  constituents.  In  particular,  composites  based  on  electroactive  polymers  have  shown  tremendous 
promise  due  to  their  actuation  capabilities  resulting  from  the  ability  to  produce  strains  of  up  to 
40%  [28,26],  and  more  than  100%  in  elastomer  actuation  [36].  The  attractiveness  of  polymers  is 
not  only  due  to  their  electromechanical  properties,  but  also  due  to  their  lightweight,  short  response 
time,  low  noise,  durability  and  high  specific  energy.  Proper  optimization  of  their  fine-scale 
structure  can  lead  to  a  vast  improvement  in  electromechanical  coupling  [26]. 

On  the  modeling  front,  there  are  numerous  nonlinear  models  [37,14,  38,39,40,41,42,43,44,45,46, 
47,48,49]  capable  of  describing  electromechanical  behavior  for  large  deformation  problems. 
There  are  a  number  of  constitutive  phenomenological  and  micromechanical  models  for 
electromechanical  materials,  such  as  ferroelectric  switching  [50,51,52,53,54,55].  So  while  the 
behavior  of  different  phases  is  well  understood  and  can  be  described  sufficiently  well,  the 
computational  cost  of  resolving  fine-scale  details  is  very  high  and  is  often  beyond  the  capacity  of 
modern  computers.  Thus  models  capable  of  describing  the  overall  behavior  of  electroactive 
composites  without  resolving  fine-scale  details  are  needed.  The  phenomenological  equations 
describing  the  behavior  of  electroactive  composites  are  usually  very  complex  and  often  unknown 
and  thus  there  is  increasing  need  for  multiscale-multiphysics  based  methods. 

There  are  a  number  of  publications  which  make  use  of  homogenization  theory  for 
electromechanical  coupling  in  composites  and  polycrystalline  materials,  but  mostly  restricted  to 
linear  piezoelectric  coupling  [1,8,56,57,58,59,60,61,62,63],  For  mechanical  problems  a  number  of 
homogenization  methods  accounting  for  large  deformations  and  arbitrary  nonlinear  constitutive 
relations  have  been  developed  [64,65,66,67,68,69,70,71,72],  In  the  present  manuscript  we  extend 
the  nonlinear  homogenization  framework  proposed  in  [65]  to  nonlinear  electromechanical 


materials.  To  the  authors’  knowledge  this  is  first  attempt  to  develop  a  general  nonlinear 
mathematical  homogenization  theory  for  electromechanical  materials.  The  derivations  of  this  paper 
are  applicable  to  arbitrary  fine-scale  structures  and  arbitrary  constitutive  models  of  phases. 


2.  Mathematical  model 

First  we  present  the  equations  describing  the  behavior  of  a  deformable -polarizable  body,  subjected 
to  mechanical  and  electrical  excitations.  The  formulation  follows  the  works  of  Yang  and  Tiersten 
[37,14].  The  approach  is  based  on  full  electromechanical  coupling  in  a  solid  via  Maxwell  stress, 
general  nonlinear  constitutive  equations,  large  deformations  and  strong  electric  fields.  Quasistatic 
approximation  is  adopted. 

Consider  deformable-polarizable  body  occupying  volume  Q  with  boundary  dfl .  When  the  body 
is  placed  in  electric  field,  differential  material  elements  experience  both  body  forces  and  couples 
due  to  electric  field.  There  is  a  full  coupling  through  the  Maxwell  stress  and  via  the  constitutive 
equations  [14,50].  Following  Tiersten  [37],  the  electric  body  force  F£ ,  couple  CE  and  power  wE 
are  used  to  derive  the  balance  equations 


Fj  =  PeEj  +  P/Ejj, 

(2.1) 

Cf  =  eijkPjEk , 

(2.2) 

wE  =pmEini, 

(2.3) 

where  Ej  denotes  electric  field,  pm  the  current  mass  density,  pe  the  current  free  charge  density 
and  7r  =  P  /  p,  the  polarization  per  unit  mass.  From  (2.1)  and  (2.2),  the  governing  equations  can 
be  derived  following  [14,37]: 

Gauss  law 


DU-Pe=  0» 


where  D  denotes  the  electric  displacement, 


(2.4) 


D.  =  s()Ei  +  P,  (2.5) 

Pt  and£'0  are  polarization  and  electric  constant,  respectively.  Since  divergence  of  polarization 
determines  induced  charge  density,  the  total  charge  density  is  given  by 


P>e=Pe-Pi,r 


(2.6) 


The  electric  field  is  vortex  free  which  can  be  written  as 


(2.7) 


£A  =  °’ 

which  implies  that  electric  field  can  be  expressed  in  terms  of  the  gradient  of  electric  potential 


The  continuity  equation  is  given  by 


Et=-*r 


Pn,  +  PmVi.i  =  0 


and  the  linear  momentum  balance  equation  is 


auj+  P, 


,b,+FiE  =  ( 


=  1  (7ii,j  +  ^j)  +  bi=Pmvn 


(2.10) 


where  cr  is  Cauchy  stress,  6  the  body  force  (excluding  electric  body  force);  erf  is  electrostatic 


stress,  which  is  related  to  electric  force  off  =  FE  . 


The  electrostatic  stress  erf  can  be  expressed  as 


erf  =  D,E,  — E  S„  =  P,E,  +  en  E,E,  — EhEh8» 


2  o  k  k  y  ~  >  J  >  J  2  k  k  ij 


(2.11) 


Divergence  of  is  given  by 


<j=DjJEi+PjEiJ=PeEi+PJEiJ=FE. 


(2.12) 


The  sum  of  Cauchy  stress  cr  and  the  electrostatic  stress  erf  gives  the  total  stress  r..  =  cr  +  erf  , 
which  is  symmetric  and  can  be  decomposed  into  symmetric  tensor 


cr  =  cr.  +PE  =  cr 

j  v  >  J  J' 


(2.13) 


and  Maxwell  stress  tensor 


of  =en  \  E,E. - E,E,S.  =cr". 


ij  \^i^j  ~^k^k”ij  ^  ji 


(2.14) 


Using  the  total  stress,  the  balance  of  linear  momentum  (2.10)  can  be  rewritten  as 


TijJ+Pn,bi=Pjr 


(2.15) 


The  angular  momentum  balance  equation  is  given  by 


Fjk^ik  +  Cf  =  si]k  (crjk  +  PjEk )  =  (o-  ,  +  ( tJ )  =  0, 


(2.16) 


or 


(2.17) 


which  suggests  that  the  total  stress  is  symmetric. 
The  conservation  of  energy  is  given  by 


p  e  =  (t  v  .  +  p  E  tc  . 

r  m  IJ  l,J  r  m  l  l 

(2-18) 

The  free  energy  can  be  introduced  through  Legendre  transform 

1 

II 

(2.19) 

and 

p  w  —  a  v  -PE  . 

r  niT  ij  i,j  i  i 

(2.20) 

The  free  energy  provides  the  constitutive  relations 

(2.21) 

\ySa). 

(2.22) 

In  the  present  manuscript  the  weak  form  of  governing  equations  will  be  expressed  in  the  initial 
configuration.  Therefore,  we  will  transform  the  strong  form  of  governing  equations  into  the  initial 
configuration.  The  balance  laws  in  the  so-called  in  Lagrangian  description  are  given  by 


(Dj,j-Pe  =  0,  (2.23) 

where  (DL  =  JFL '  D.  and  pE-peJ  are  Lagrangian  description  of  electric  displacement  and  charge 
density,  Z>  =  L<DL . 

Electric  field  is  vortex  free  in  the  Lagrangian  description  and  is  given  by 

£,A,j  =  0,  (2.24) 

where  (EK  =  EiFiK  =  ~<f)AFi  K  =  -<j)K  is  Lagrangian  description  of  electric  field,  Ei  =  (EK  FK\ . 
Equation  (2.24)  can  be  used  implicitly  to  describe  electric  field  as  a  gradient  of  electric  potential. 

KkL,L+Bk=  PMvk,  (2-25) 

where  Ku  =  JFL  '  r(/  and  pM  -  pmJ  are  the  first  Piola-Kirchhoff  stress  and  material  density, 
respectively,  and  r;/  =  J  ]Fj  ,KU  .  KiL  is  nonsymmetric  satisfying 


£kijFj'LKiL  =  0, 


(2.26) 


which  confirms  that  moments  resulting  from  electric  body  force  acting  on  the  infinitesimal  volume 
are  equilibrated  by  internal  forces. 

The  energy  balance  equation  is  given  by 


PmV  =  TsklSkl-Vk‘Ek,  (2.27) 

where  the  second  Piola-Kirchhoff  stress  T^L  is 

T^=JFK1kFL~ykl.  (2.28) 

The  Green-Lagrange  strain  is  given  by 

^KL  =  ( UK,L  +  UL,K  +  UM,KUM,L  )  ^  (2.29) 

and  the  Lagrangian  description  of  polarization  vector  <PK  is 


VK  =  JFkA-  Pt=J% KCPK,  (2.30) 

^=ai/  +  ^j=J  'FlkFu7%.-  (2.31) 


The  energy  balance  determines  constitutive  equations,  which  can  be  written  as 

dy/ 


tskl=t^l(skl,<ek)  =  p, 

F  K  =  F  K  (  ^kl  ’  F K  )  =  ~Pm 


dSa 

dy/ 

dEr  ’ 


Kn  and  <Dl  may  be  calculated  from 


dy/ 


K,l=FlkPm- - +  JF,S, 


dSk 


1°  0 


EE.--E.E.S 

i  j  ^  k  k  ij 


®  k  £oFFklFl  Pm 


dy/ 

d<E. 


(2.32) 


(2.33) 


To  complete  the  strong  form  (2.23),  (2.25),  (2.32),  (2.33)  of  the  initial-boundary  value  problem  is 
given  by 


ui  =  w(  on  dCl“ , 

(p  =  (p  on  dO? , 
KkLNL=TK  on  5Qr, 
DkNk  =  -do  on  dQ? , 
n(  =  ui  (0)  at  t  =  0, 


duj 

dt 


(2.34) 


where  co  is  surface  charge  density  applied  on  the  portion  of  the  surface  SQ*  and  TK  is  mechanical 
traction  applied  on  the  portion  of  the  surface  5Q7  ;  ui  and  <p  are  displacements  and  electric 
potential  prescribed  on  the  surfaces  <9Q"  and  dO.<p ,  respectively,  such  that 

8Q.W  n  dQ.*  =  dO!‘  n  =  5Q 

dQTvjdQ?  =aQ"u5Qr  =0 


(2.35) 


The  free  energy,  which  determines  constitutive  relationships  for  nonlinear  electroelastic  materials, 
can  be  written  as  power  series  of  SAB  and  (EA  [14] 

PmV^KL^k)  ABCD^AB^SD  ~eABCC^’A^BC  ~  ^  X  2  AB^A^B  +  ^  C2  ABCDEF^  AB^ CD^ EF  + 


+  2  ^ABCDE^A^ BC^DE  ^  ^ ABCD^A^B^ CD  r  XlABC^A^B^C  +  ^  „  C4ABCDEFGH^ AB^CD^ EF^GH  + 


1 

2 


1 

6' 


1 

24 


(2.36) 


+  -k. 


1 


1 


1 


2ABCDEFG 


^■A^BC^DE^FG  a\ ABCDEF^A^B^CD^EF  s-  ^iABCDE^A^B^C^ ' DE  .  %4 ABCD^A^B^C^, 


where  constants  c1ABCD ,  eABC ,  Ziab’  c2abcdef  >  k abcde  >  P abcd  >  Xaabc  >  c4 abcdefgh  >  ,v2  abcdefg  > 


,  A 


IABCDEF 


,  kiABCDE  X 4 abcd  are  often  referred  to  as  the  fundamental  material  constants.  The  structure 


of  y/^S^,1 Ek)  depends  on  the  symmetries  of  particular  materials  under  consideration  [73,74]. 


From  the  above  strong  form  (2.34)  it  follows  that  the  electric  fields,  <£  and  <P  ,  are  coupled  to  the 
mechanical  fields,  F  and  K  ,  at  two  levels  [50],  First,  is  the  electrostatic  coupling,  when  the 
electric  fields  directly  generate  distributed  forces  via  Maxwell  stress,  which  in  turn  affects  the 
mechanical  equilibrium  of  a  solid.  The  resulting  stresses  depend  on  the  second  order  terms  of 
electric  field,  and  while  typically  very  small  (1(T5  MPa  for  a  lMV/m  field),  could  have  a  significant 
effect  if  the  fluctuations  of  the  electric  field  are  large  (field  singularities)  as  in  the  case  of 
heterogeneous  media  with  high  contrast  in  the  dielectric  moduli  or  in  electrodes  and  conducting 
crack  tips  [75,76,26],  This  nonlinear  phenomena  is  universal  and  common  for  both,  dielectrics  and 
conductors  [77],  and  has  no  converse  effects,  since  the  mechanical  stress  state  does  not  directly 
influence  the  electrostatic  balance  [50],  but  can  affect  it  through  the  motion,  when  such  a  motion 
changes  electric  field. 

The  second  level  of  coupling  occurs  in  the  constitutive  behavior  of  dielectric  materials,  expressed 
by  equations  (2.21),  (2.22)  or  (2.32),  (2.33).  Polarization  induces  strain,  while  mechanical  stress 
affects  polarization,  which  introduces  full  coupling.  Most  common  electromechanical  materials  are 
piezoelectrics,  ferroelectrics  and  electrostrictors  [78], 

In  piezoelectrics,  there  is  a  linear  dependence  between  applied  stress  and  induced  electric 
displacement  (direct  piezoelectric  effect) 


A  =  dm(TJk 


(2.37) 


and  between  applied  electric  field  and  induced  strain  (inverse  piezoelectric  effect), 


where  djk  are  piezoelectric  coefficients.  A  positive  electric  field  results  in  positive  strain,  and  vice 
versa,  negative  electric  field  results  in  negative  strain. 

Ferroelectric  materials  have  spontaneous  polarization  Ps  in  absence  of  external  electric  field.  The 
charge  due  to  spontaneous  polarization  is  usually  masked  by  the  charges  from  surroundings  of  the 
material.  Direction  of  spontaneous  polarization  can  be  switched  by  external  electric  field.  Thus 
ferroelectric  polycrystallines  are  characterized  by  hysteresis  in  polarization-electric  field  P(E)  and 
induced  strain-electric  field  e(E). 

In  electrostrictive  materials  induced  strain  is  proportional  to  the  square  of  polarization 

£ij  =  Qykipkpn  (2-39) 

where  Qyki  are  electrostrictive  coefficients.  For  mildly  strong  electric  fields,  the  dependence 
between  polarization  and  electric  field  is  linear  and  (2.39)  may  be  written  as 

£<  =  MjkIEkE„  (2.40) 

where  Mijki=ZkmZinQkmin  and  Zij  is  dielectric  susceptibility.  Formulation  of  ‘converse’ 
electrostrictive  effect  is  also  possible  [79].  Here  both,  positive  and  negative  electric  fields  result  in 
positive  longitudinal  strain. 

In  both,  electrostriction  and  electrostatic  coupling  the  overall  behavior  is  detennined  by  the 
average  of  the  fluctuation  of  the  electric  field  (rather  than  its  average).  Therefore  one  can  obtain 
large  coupling  with  relatively  small  macroscopic  field  [19,35].  Huang  et  al.  [35]  described  a  three- 
phase  polymer  based  actuator  with  more  than  8%  actuation  strain  attained  with  an  activation  field 
of  20  MV/m. 


3.  Mathematical  homogenization  for  nonlinear  coupled  electromechanical  problem 

In  this  section  we  will  introduce  scale  separation  and  proceed  with  deriving  equations  for  different 
scales  from  the  strong  governing  equations  in  the  Lagrangian  description  derived  in  the  previous 
section. 

Consider  a  body  with  volume  CiA  and  surface  cTV  made  of  heterogeneous  solid  with 
characteristic  size  of  the  heterogeneity  /  which  is  much  smaller  than  the  body  dimension  so  that 
heterogeneities  are  considered  to  be  “invisible”.  The  body  is  effectively  homogeneous  at  the  coarse 
scale;  its  fine-scale  structure  can  be  observed  by  zooming  (Figure  1)  at  any  coarse-scale  point  A. 
We  will  denote  the  volume  and  boundary  of  the  body  on  the  coarse-scale  by  Q  and  dQ , 
respectively. 


To  describe  the  dependence  of  various  fields  VP/  (X)  on  the  coarse-scale  and  fine-scale 

coordinates  we  will  introduce  global  coordinate  system  X  and  local  (unit  cell)  coordinate  system  Y 
associated  with  fine-scale  structure  placed  at  every  coarse-scale  point.  The  two  coordinates  are 
related  by  Y  =  X//l,  0<T<scl.  Dependence  of  the  field  X¥A  (X)  on  the  two  scale  coordinates  is 
denoted  as 


'Pi(X)  =  'T(X,Y)  (3.1) 

where  superscript  X  denotes  existence  fine-scale  details. 


We  assume  that  the  body  is  composed  of  a  periodic  repetition  of  unit  cells  (UC)  with  volume  0y . 

The  periodicity  may  be  global,  when  the  whole  body  is  a  lattice  consisting  of  unit  cells  0K  (Figure 

2. a),  or  local,  when  periodicity  holds  only  in  a  small  neighborhood  of  coarse-scale  point  (Figure 
2.b).  Unit  cell  coordinates  Y  in  (3.1)  are  defined  with  respect  to  the  initial  (undeformed)  fine-scale 
configuration. 


d, Q;- 

(a) 


(b) 


Figure  2.  Global  (a)  and  local  (b)  periodicity. 


To  construct  the  weak  form  of  the  governing  equations  we  will  need  to  integrate  over  the  volume 
i¥  .  This  can  be  carried  out  as  long  as  the  size  of  unit  cell  is  infinitesimally  small,  which  yields  the 
following  fundamental  lemma  of  homogenization  ([80,81] 


Jim  1  'r(X,Y)rfn=  limjU|  J  V(X,Y)d® 

nl  nUwr|©T 


d  Q, 


(3.2) 


which  also  implies  that  the  fine-scale  domain  0y  exists  at  every  point  in  the  coarse-scale  domain 
(Figure  3.). 


Figure  3.  Transition  from  the  heterogeneous  to  homogeneous  domain  as  A  — >  0+ 

To  derive  the  coarse-scale  equations  we  will  start  from  the  equilibrium  equation  and  Gauss  law  in 
the  Lagrangian  description  (2.23),  (2.25).  As  before  superscript  A  denotes  dependence  of  a 
quantity  on  fine-scale  details 

K!L,L(F\-EA)+Bi(X)=piv,  (3.3) 

<D‘j(f’-^)-p1e(X)  =  0  (3.4) 

with  initial  and  boundary  conditions 


uj  =  ui  on  3Q" , 

(p  =  (f>  on  dO!p , 
KkLNL=TK  on  dQr, 
DkNk  =-(o  on  3Q®, 
uf  =  ui  (0)  at  t  =  0, 

=  0)  at  t  =  0, 


(3.5) 


dClAw  =  dQ.Au  n5Qir  =  SO, 

=  dQAtl  udQAT  =  0 


(3.6) 


Utilizing  the  asymptotic  expansion  [82,83,84]  of  displacement  and  electric  potential  fields  yields 

uA  (X)  -  u,  (X,  Y)  =  u°  (X)  +  Au\  (X,  Y)  +  Z2u:  (X,  Y)  +  o(A3 ) ,  (3.7) 

cpx  (X)  =  <p(X,  Y)  =  cp°  (X)  +  V  (X,  Y)  +  A2(p2  (X,  Y)  +  0(; l3 ).  (3.8) 

The  size  of  the  unit  cell  is  of  the  order  A  and  in  the  limit  it  is  assumed  to  be  infinitesimally  small, 

so  that  first  terms  in  asymptotic  expansions  for  uf(X)  and  <p"  ( X )  do  not  depend  on  fine-scale 

coordinate  Y.  This  has  been  shown  to  be  a  unique  solution  for  linear  elliptic  problems.  Note  that 
for  certain  class  of  nonlinear  problems,  such  as  problems  involving  softening  and  localizations  as 
well  as  neutronic  diffusion,  radiative  transport  problems  [85],  one  has  to  consider  large  oscillations 
in  the  leading  order  term  in  which  case  the  first  term  in  asymptotic  expansion  will  depend  on  fine- 
scale  coordinates  Y, 


Expanding  every  term  in  (3.7)  and  (3.8)  in  Taylor  series  around  the  unit  cell  centroid  X  =  X  yields 

u°  (X)  =  u°  (x\  +  A^—  Y;+y i2  8  U<  Y,Yk+0(A 3),  (3.9) 

lK  ’  '  1  >  dXj  .  J  a XjdXK  ±  K  ’ 

n  n 

u"  ( X,  Y )  =  n”  (x,  y)  +  T  —X  Y.+A2  1  Y.Yk+0(A3),  (3.10) 

'  V  1  >  dXJx  dXjdXK  ±JK  1  ’ 

/(X)  =  p»(x)+A +  F,F,  +o(X=),  (3.11) 

«>‘(X,Y)  =  «.-(X,Y)  +  A|^  Y,+A2- AX—  F,FA.+o(l5),  (3.12) 

OAj  k  OAjCAk  k 

where  Xj-Xj  =  AYj  was  used. 

Substituting  the  above  into  asymptotic  expansions  (3.7),  (3.8)  gives  the  modified  asymptotic 
expansion 

uA  (X)  =  ut  (X,  Y)  =  u°  (x)  +  Au)  (X,  Y)  +  A2u:  (x,  y)  +  O (A3 ) ,  (3.13) 


where 


;,"(x)=«»(x),  «;(x,y)=«;(x,y)+ 


du, c 


dX, 


(x,  y)  =  u]  (x,  y)  + 


3m,1 


dX, 


2..0 


Yj+—  ^ 

x  2  dXjdXK 


Y,, 


Y  Y 

1  J1  k 


(3.14) 


and  similarly  for  potential 


<P‘ 


(X)  =  (p(x,  Y)  =  ft  (x)  +  A  ft  (x,  y)  +  A28p 2  (x,  y)  +  0(A3),  (3.15) 


where 


<P 


<P' 


(x)  =  ft (x),  ft (x, y)  =  ft (x, y)  + 


8  ft 


8X , 


Y, 


(x,y)  =  ^(x,y)+ 


dft 


ex, 


1  82ft 


Yj  + 

x  2  dXjdXk 


YjYk 


(3.16) 


dfA  df(X  Y)  1  8f(X  Y) 

In  the  classical  homogenization  the  spatial  derivative  is  defined  as  — —  =  - — - — - — -h - : — - — - — 

8Xf  8XI  c  8YI 

.  Since  in  the  modified  asymptotic  expansions  (3.13)  and  (3.15)  all  tenns  are  calculated  with 
respect  to  the  UC  centroid  X  ,  there  is  no  dependence  on  X  and  the  spatial  derivative  reduces  to 


dfl(X)  _  1  5/(x,y) 

8Xt  ~  A  3Y(. 

Thus  displacement  gradient  may  be  expressed  as 

l  3m,(x,y)  3w!(x,y)  3m2(x,y) 

—  -  +  A- 


(3.17) 


3m; 


dXK  A  8Yk 


37, 


8Y, 


o(a2), 


(3.18) 


where 


3m;(x,y)_3m;(x,y)  ^  gM°(x) 


37, 


37, 


dXt 


3m;(x,y)_3m2(x,y)  5m;(X,Y) 
8Yk  ~  8Yk  8Xk 


82u](X,  Y) 


.  +  8X,8YK 

X  J  K 


Yj  + 


32m;°  (X,  Y) 


(3.19) 


8Xj8Xk 


7. 


The  asymptotic  expansion  of  the  deformation  gradient  is  given  as 


fX‘>iK+X-=fiI(x,y)+zf,'k(x,y)+°(z2), 


(3.20) 


where 


»;(x,Y)  0u,‘(x,y)  a„»(x) 

Fi(x,Y)  =  Sx+— =  <?„-+  '  •'  > 


31' 


31' 


-  +  - 


dXr 


=  ^(X,Y)+^(X,Y),  (321) 


F; 


iK 


(*H 


-  + 


3m, 0  (X) 


3XC, 


3x, 


3*V 


>  FiK 


(X,Y): 


a'.‘(*>Y)  irl^vU  asfM 


3Y, 


4(x,y) 


37 


Similarly,  the  electric  field  is  given  as  negative  derivatives  of  electric  potential 


where 


jy  ’  dXj  X  dYj  8Yj  dYj  1  ’ 


=  <Ej  (x,  Y)  +  A<e]  (x,  y)  +  0(A2 ) , 


(3.22) 


(3.23) 


<E 


,°(X,Y): 


8q>1  _  V(X,  Y)  3^°(X) 


37, 


37, 


3X, 


:<(x)  +  e;(x,y),  (3.24) 


E/(X) 


V(x) 


3X, 


, *;(x.y)^M  e;(x,y)=  ^ 


37, 


(3.25) 


The  coarse-scale  deformation  gradient  f£  and  electric  field  <E</  are  related  to  the  average  of  the 
leading  order  defonnation  gradient  F°K  (x,y)  in  (3.21)  and  electric  field  <E°  (x,  y)  in  (3.24)  over 
the  unit  cell  domain 


„  /  -  \  l  p  o/~  \  3m, 0  ( X) 

^  (x)  =  ^7  J  (x,  Y)rf0,  =  4,  + 

Y  @K  K 


(x)  =  —7  J  ‘Ej  (x,  y) 

I  @r 


dcp° 


dx, 


(3.26) 

(3.27) 


c)u^ 

where  conditions  [  — - d 0  =  0  and  [  <^— d 0  =  0  were  employed.  These  conditions  correspond 

J  p)V  J  r)Y 


d(f »' 


37 
e,  u1j 


37 


to  periodicity  (or  weak  periodicity)  of  fine-scale  displacement  and  potential  fields,  which  will  be 
discussed  in  the  next  section. 


Expanding  stress  and  electric  displacement  in  Taylor  series  around  the  leading  order  deformation 
gradient  F°(X,Y)  and  electric  field  (£°(x,y)  yields 


K*(f\<Ez)  =  Ku{f°  ,(E°)  +  2- 


FL+^ 


■El+0{r‘)  =  K;^ZK‘J+0(22) 


mN  |p°  £° 


“’a/  If0  e° 


(3.28) 


( D1(f\<ea)  =  (Dj(f0,‘e°)  +  a - 


'J  Z7l 


mN  Ip0  e° 


pi  |  ;  d(DJ 

M  SE„ 


£^+o(A2)  =  ©>/1©)+o(/12)  (3.29) 


’M  Ip0  e° 


Taylor  expansion  of  stress  and  electric  displacement  around  the  unit  cell  centroid  X  yields 

X„(f\£j)  =  X„(X,Y)  =  X,;(x,y)  +  a[|S  Y.+Kl (x,y)|  +  0(22),  (3.30) 


©,(fj,^)  =  ©j(x,y)=©;(x,y)  +  1  — i  Yk+<D){%  Y)  +0(22).  (3.31) 

yoxK  . 

Again,  here  the  relation  (x,.  -Xj  j  =  AYt  has  been  utilized.  Further  substituting  above  into 
equilibrium  equation  (3.3)  and  Gauss  law  (3.4)  yields 


fflf,«(x,Y)  *01  a^(x,Y) 


+s,(x,y)-p„(x,y) 


d2w,°(x,t) 


+e3m_pt%  Y)=0. 

r)Y  r)Y  r)Y  E  \  > 


(3.32) 


(3.33) 


Collecting  terms  of  equal  orders  in  (3.32),  (3.33)  yields  the  two-scale  equilibrium  equation  and 
Gauss  law 


o(r'): 


ax“(x,  y)  s©°(x,y) 


(3.34) 


0(X°): 


Wj  dYJ 

dV]  5©j(x,Y) 


5,(x,y)  =  A/(x,y) 
=  ^(x,y). 


82u°  (x,t) 


(3.35) 


Integrating  0(At>]j  equations  over  the  unit  cell  domain  and  exploiting  J 

© 

0©j(x,Y) 


SK'U  (x,Y) 


dY, 


dY  =  0 . 


dY, 


-dY  =  0 ,  which  correspond  to  periodicity  (or  weak  periodicity)  of  fine-scale  stress  and 


electric  displacement,  yields  the  coarse-scale  equations 

a*£(x)| 


ax, 


f2n,°(x,t) 

+Bc\(i 

X 

II 

<  \  ’ 

1  dr 

8©t(X) 


dX, 


(3.36) 


=  p£c(x), 


where 


X-u  (x)  =  1  X^(F°(x,  y),e(x,  Y))^@,  *,c(x)  =  T-  Jii,(x,Y)rf©,  (3.37) 

Qy  Qy 

©,c(x)  =  -|-|©;(fo(x,y),£(x,y))J0,  pcE  (X)  =  — j—  J  pE  (x,  Y^©.  (3.38) 

By  By 

Hereafter  we  will  make  use  of  an  argument  of  the  unit  cell  infinitesimality  by  which  the  unit  cell 
centroid  X  can  be  positioned  at  an  arbitrary  point  X .  Therefore,  the  centroid  X  can  be  replaced 
by  an  arbitrary  point  X  .  Equations  (3.36)  with  boundary  conditions  (3.5)  define  the  coarse-scale 
boundary  value  problem. 

The  set  of  equations  (3.34)  together  with  periodic  (or  weakly  periodic)  boundary  conditions, 
discussed  in  the  next  section,  defines  fine-scale  boundary  value  problem. 


4.  Boundary  conditions  for  the  unit  cell  problem 

In  this  section  various  boundary  conditions  imposed  on  the  unit  cell  will  be  considered.  The  most 
common  are  the  periodic  boundary  conditions,  when  the  displacement  and  electric  potential 
perturbations  on  the  opposite  sides  of  the  unit  cell  are  assumed  to  be  equal. 

Consider  the  asymptotic  expansions  for  displacement  and  electric  potential  fields  (3.13),  (3.15) 


{Xd)  =  u\%Yp)  =  u^[%t)+AUF^[±d)-Suyj+u)[±,Y,t)\  +  0{X1), 


^(X,i)  =  ^(x,Y,i)  =  ^°(x,i)  +  ^r-iEj(x,i)7y+^1(x,Y,i) 


+o(a2), 


where  (3.25)  and  (3.26)  were  used.  The  leading  order  terms,  which  represent  the  rigid  body 
translation  of  the  unit  cell  and  constant  potential,  are  independent  of  the  unit  cell  coordinates  Y. 
The  O  ( A )  terms  describe  the  unit  cell  distortion  and  electric  field  in  it.  The  terms 

^  /yj  |  X,  t  j  -  j  and  <Ej  (  x,  /  j  T;  represent  a  uniform  coarse-scale  deformation  and  uniform 

coarse-scale  electric  field,  while  the  terms  and  (f  (  x,  Y,/  j  capture  the  deviations  from 

the  uniform  fields  induced  by  material  heterogeneity.  These  terms  are  important  for  electrostriction 
effect  [26].  Making  these  terms  larger  increases  the  electrostriction  effect. 

Figures  4(a)  and  4(b)  show  the  initial  and  defonned  shape  of  the  unit  cell,  respectively.  The  dotted 
line  in  Figure  3(b)  depicts  the  deformed  shape  of  the  unit  cell  due  to | F.j  (X,/j  -d';/  j  Yj ,  whereas 
the  solid  line  shows  the  contribution  of  the  two  0(A)  terms. 

At  the  unit  cell  vertices  dQfi  ,  the  deviations  from  the  unifonn  fields,  u)  (x,  Y,/  j  and  cp'  (x,  Y,/ j , 
are  assumed  to  vanish.  For  the  remaining  points  on  the  boundary  of  the  unit  cell,  the  deviations 
from  the  average  u)  (X,  Y,/ j,  <p)  ^X,  Y,/j  are  prescribed  to  be  periodic  functions.  Figure  4  depicts 

two  points  M  and  S  on  the  opposite  faces  of  the  unit  cell  with  M  and  S  tenned  as  the  master  and 
slave  points,  respectively.  The  displacements  of  the  two  points  are  given  as 

«;(x,Y  “  ,t)  =  [Ff%t)-Suyy  +«'(x,Y (4.1) 

fi,1(x,Ys,()  =  (^(x,()-4)l'/+U;(x,Ys,<).  (4.2) 

Subtracting  equation  (4.2)  from  (4.1)  and  accounting  for  periodicity,  i.e. 
u)  ( X,  Y  w ,  /  j  =  u)  (X,  Ys,/| ,  gives  a  multi-point  constraint  (MPC)  equation 

«‘(x,Y",()-.:‘(x,YI,()=(^(x,/)-4)(if  -rf)  (4.3) 

where  FM  and  Ys  represent  the  local  coordinates  of  the  master  and  slave  nodes  on  the  unit  cell 
boundary,  respectively. 

Similarly,  periodicity  of  electric  potential  perturbations  cpx  fx, \M ,t\  =  (pl  (x,  Ys,tj  gives  the 
following  constraint  equation 

v{X,Yu ,t)-p%Ys ,t)  =  -y(x,t)iyY;<  -Yp. 


(4.4) 


(a)  Initial  unit  cell 
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(b)  Deformed  unit  cell 


Figure  4:  Definition  of  periodic  boundary  conditions 


Periodic  boundary  conditions  can  be  used  for  periodic  heterogeneous  medium  and  when  the  unit 
cell  distortion  is  not  considerable,  otherwise  they  are  not  physical.  For  nonperiodic  medium  or  in 
case  of  large  unit  cell  distortions  more  general  boundary  conditions  have  to  be  used  instead.  Note 
that  the  only  time  when  the  periodicity  condition  was  exercised  was  in  deriving  equations  (3.26), 
(3.27).  For  (3.26),  (3.27)  to  hold  the  following  conditions  must  be  satisfied 


f  r  =  0,  r 

}@r  qYj  }&y  dYj 


(4.5) 


Applying  Green’s  theorem  and  exploiting  relations  (3.14)b,  (3.16)b,  and  (3.26),  (3.27)  the  above 
reduces  to 


1  (",l(x,Y,/)-(^(x,()-4.)FA.)^.rfr,  =  o, 
|  (x,  Y,  t )  -  (-<  (x,  t ))  Yk  )  NKdTY  =  0, 


where  50  y  is  the  boundary  of  0y  and  N  K  are  the  components  of  the  unit  normal  to  the  boundary 
50y .  Equation  (4.6)  represents  the  so-called  weak  periodicity  condition. 


Alternatively  to  equations  (4.3),  (4.4)  and  (4.6)  an  essential  boundary  condition 

^(x,Y,f)-(-<(x,f))rA.  =o 


(4.7) 


is  often  exercised  in  practice.  It  corresponds  to  zero  perturbations  from  coarse-scale  fields  on  the 
unit  cell  boundary. 


The  essential  boundary  conditions  can  be  enforced  either  in  the  strong  form  (4.7)  or  in  the  weak 
form  as 
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where  Juj ,  /uip  are  the  Lagrange  multipliers  representing  unknown  tractions  and  surface  charge 
density  on  d&Y.  If  we  choose  //  =  K(uNj,  nv  =  where  K(u ,  (D^  are  constant  over  0}  and 

require  (4.8)  to  be  satisfied  for  arbitrary^  and  ,  we  then  obtain  equation  (4.6).  Therefore, 

equation  (4.6)  will  be  referred  to  as  a  weak  compatibility  boundary  condition  while  equation  (4.7) 
as  the  strong  compatibility  condition.  Equation  (4.6)  is  in  the  spirit  of  the  work  of  Mesarovic  [86] 
who  imposed  the  unit  cell  to  satisfy  average  small  strains. 

The  different  boundary  conditions  can  be  denoted  for  generality 
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=  0  on  50. 
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(4.9) 


5.  Two-scale  problem 

The  two-scale  problem,  consisting  of  the  unit  cell  equations  (3.34)  subjected  to  periodic  (or  other) 
boundary  conditions  (4.9)  and  the  coarse-scale  equations  (3.36),  is  two-way  coupled.  The  link 
between  the  two  scales  is  schematically  shown  in  Figure  5. 


coarse- scale 


D  VP 


Figure  5:  Information  transfer  between  the  coarse-scale  and  fine-scale  problems 

The  fine-scale  problem  is  driven  by  the  overall  (coarse-scale)  deformation  gradient  Ff  (  x,/ j  and 

electric  field  F(K  |  X,  /  j .  They  are  calculated  for  every  material  point  of  the  coarse-scale  problem 

(in  practice  only  for  integration  points  of  the  coarse  mesh)  and  used  together  with  (4.9)  to 
formulate  boundary  conditions  to  be  imposed  on  the  unit  cell.  Once  the  unit  cell  problem  is 
solved,  it  provides  the  coarse-scale  problem  with  coarse-scale  stress  Kcu  and  electric  displacement 

Qj  via  (3.37),  (3.38). 

This  effectively  provides  a  coarse-scale  constitutive  relationship.  Additionally,  the  local  coarse- 
scale  consistent  tangent  is  derived  from  the  fine-scale  stiffness. 

This  scale-bridging  approach  belongs  to  the  category  of  information-passing  (sometimes  referred 
to  as  hierarchical  or  sequential)  multiscale  methods  which  evolve  a  coarse-scale  model  by 
advancing  a  sequence  of  fine-scale  models  in  small  windows  (representative  volume  or  unit  cell) 
placed  at  the  Gauss  points  of  the  discretized  coarse-scale  model. 

The  coupled  two-scale  problem  is  summarized  below: 


(a)  Fine  scale  problem 


(5.1) 


(b)  Coarse-scale  problem 


(5.2) 


6.  Finite  element  discretization 

Both  the  coarse-  and  fine-scale  problems  are  discretized  using  finite  elements.  The  displacement 
and  electric  potential  fields  of  the  fine-scale  problem  u)  (x,  Y,/ j ,  ip'  (X,  Y,/ j  are  approximated  as 

u;(x,Y,tWl(Y)4(xA 

;  ;  (6.1) 

^(x,Y,t)  =  4(Y)4(x,t), 

where  subscript  B  denotes  the  node  number,  4(F)  are  the  unit  cell  shape  functions  and  d)B,  <p\ 
are  the  nodal  displacements  and  potentials  in  the  unit  cell  mesh.  Let  d™  (x,f),  4  (x,t)  be  the 


master  (independent)  degrees  of  freedom  and  express  d\R ,  <p\  by  a  linear  combination  of  dR.  ( X,/ 
and  tp™  (x,  /  )  defined  by 
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so  that  the  constraint  equation  (4.9)  in  the  discrete  form  is  satisfied 
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(6.3) 


Then  writing  the  Galerkin  weak  form  of  (5.1)  and  discretizing  it  using  (6.1)  yields  the  discrete 
residual  equation: 
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where  the  left  subscript  and  superscript  denote  the  load  increment  and  the  iteration  count  (for 
implicit  method)  at  the  coarse-scale,  respectively;  H\rfR  and  r* 1 ,  Ad'R  and  Ai/f  are  the 

residuals,  displacement  and  potential  increments  in  the  (m  +  l)lh  iteration  of  the  (n  +  l)th  load 
increment,  respectively.  If  the  constitutive  equations  are  defined  in  terms  of  the  Cauchy  stress  and 
electric  displacement  in  the  current  configuration  it  is  convenient  to  restate  the  unit  cell  problem  as 
follows: 


Given  and  h^,„A  find  ^A^  such  that: 
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where  we  have  exploited  the  relation  between  the  first  Piola-Kirchhoff  stress  KiK ,  electric 
displacement  d)K  in  the  Lagrangian  description  and  their  current  configuration  counterparts, 
Cauchy  stress  cr.  and  electric  displacement  Di 


(6.6) 


and/ in  (6.6)  is  the  determinant  of  F jK . 


J°ij  =FjkK;k, 

JDt=FjDK 


Similarly,  the  coarse-scale  displacements  and  potential  U-(X,t),  <p°(X,t )  are  discretized  as 


u^X,t)  =  NcA{X)dcu{t), 

<p°(X,t)  =  NcA(X)<pcA(t), 

where  NA(X),  dfA  and  cpcA  are  the  coarse-scale  shape  functions,  nodal  displacements  and 
potential,  respectively.  Writing  the  weak  form  of  (5.2)  and  using  discretization  (6.7)  the  discrete 
coarse-scale  equations  can  be  expressed  as 

Given  „+15,c,,1+1pec,  nJt  and  n+1®,  find  Ii+1A d°A  and  /1+]  A^’such  that: 
n+irUC (n+l^dC ,  n+lA<pC)  =  MijAB  nJ%  +  n+j;A  -  n+J-‘  =  0 

n+sf(„+Adc,  n+Avc)  -  n+jr  -  n+jr = o  (6  8) 

n+ldU  =  n+XFU  0tl  5G" 

n+dP(A  =  „+lPA  011 

n  <—  n  + 1,  Go  to  the  next  load  increment 

where  n+lrjAc ,  n+l r'Ac  and  n+1A dfA,  n+1A (pcA  are  the  coarse-scale  residuals,  displacement  and  potential 
increments  in  the  (n  +  l)th  load  increment,  respectively,  and 
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+  r°didr, 

(6.13) 

where  MjJAB  is  the  mass,  and  the  internal  and  external  forces,  respectively.  It  is  again 
convenient  to  express  the  coarse-scale  governing  equations  in  the  deformed  configuration 


M„ab  =  S„\a  pcmNcANc,dn„  (6.14) 

f!“=l^dn-  (6J5) 

f,T  =  J„  Kbfdn,  +  Ncj,d  r„  (6.16) 

=  (6J7) 

/r  =  f  + jm.  Nfadr, ,  (6. 1 8) 

where 

pc=pcM/Jc,  bf=Bf/Jc,  pce=pcEUc,  TidTx  =  TtdT,  wxdTx  =  MT,  (6.19) 

tf=F^/Jc  (6.20) 


and  Jc  is  the  determinant  of  F(u  ;  Q.x ,  dQ.x,dQ.x,dTx  denote  volume,  infinitesimal  volume, 
boundary  and  surface  element  in  the  current  configuration. 

We  now  focus  on  deriving  a  closed  form  expression  for  the  overall  Cauchy  stress  cr f  and  electric 
displacement Zf  .  Substituting  (6.6)  into  (3.37),  (3.38)  and  recalling  c/0  =  JdQy ,  we  have 


*£(*,0  =  1^  FkIct,J®„ 
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(6.21) 


Inserting  (6.21)  into  (6.20)  and  denoting  the  volume  of  the  coarse-scale  (intermediate) 
configuration  as  I©*  I  =  J  |@y| ,  the  overall  quantities  can  be  expressed  as 
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where  AF. 1  =  F(jK FK'n  maps  the  fine-scale  deformed  configuration  0rinto  the  coarse-scale 
deformed  configuration  0*  as  illustrated  in  Figure  6. 


Figure  6:  Unit  cell  configurations:  (a)  initial,  (b)  coarse-scale  (intermediate),  (c)  fine-scale  deformed 
(final) 

The  coarse-scale  problem  may  be  solved  using  either  explicit  or  implicit  time  integration  [87]. 

7.  Numerical  examples 

In  this  section  we  will  consider  a  numerical  example  illustrating  the  ability  of  the  mathematical 
homogenization  to  resolve  the  coarse-scale  behavior  of  heterogeneous  materials  with  nonlinear 
electromechanical  coupling  subjected  to  electric  field.  The  results  of  the  mathematical 
homogenization  (to  be  referred  as  MH)  will  be  compared  to  the  direct  numeric  simulation  (DNS) 


where  a  very  fine  mesh  is  employed  to  resolve  fine-scale  details.  For  simplicity,  we  will  consider 
two-dimensional  problem  (plane  strain). 


For  nonlinear  electromechanical  material  we  will  use  a  simple  relaxor  ferroelectric  material  model 
proposed  in  [50],  where  polarization  and  strain  are  used  as  independent  variables.  It  is  assumed 
that  material  is  isotropic,  where  the  stress  depends  linearly  on  total  strain.  The  polarization  induced 
strain  sE  depends  on  the  square  of  polarization,  i.e.,  it  is  electrostrictive  material,  where 

polarizations  saturates  to  Ps  at  high  electric  fields.  The  internal  energy  function  for  relaxor 
ferroelectric  was  proposed  by  Horn  and  Shankar  as 
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(7.1) 


eE  is  polarization-induced  strain  defined  as 


^  =  Q  :  (PP)  =  ( Q ,  -  Qu )  PP  +  Q12  |P|2 1,  (7.2) 

C  and  Q  are  an  isotropic  elastic  stiffness  matrix  and  an  isotropic  electrostrictive  strain  coefficient 
matrix;  k  is  material  constant;  coefficients  Qn  and  Qn  are  defined  so  that  the  longitudinal  and 

transverse  induced  strains  relative  to  polarization  direction  are  Qn  P  and  Qn  P  “ ,  respectively. 


Stress  and  electric  field  are  calculated  from  the  internal  energy  by  differentiation 
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(7.4) 


The  first  term  in  (7.4)  represents  the  converse  electrostrictive  effect,  while  the  second  term 
represents  the  stress-free  dielectric  behavior.  Polarization  can  be  written  in  more  compact  form  as 


P  =  P5tanh(£|R|)-^-, 

rv 


(7.5) 


where 


R  =  E  +  2(f-Q:(PP)):C:Q  P  =  E  +  2<r:Q  P  (7.6) 

For  a  given  strain  and  electric  field  (which  are  known  from  the  previous  iteration),  one  can  solve 
(7.5)  for  induced  polarization  and  then  use  (7.3)  to  calculate  stress.  The  Jacobian  for  the 
constitutive  model  can  be  derived  from  (7.5)  and  (7.3)  as 
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H  =  tanh(A:|R|) 
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Y(GPa) 

V 

Qn(m4/C2) 

Qi2(m4/C2) 

Ps(C/m2) 

k(m/MV) 

115 

0.26 

1 .33x10'“ 

-6.06x10'3 

0.2589 

1.16 

Table  1.  Model  parameters  for  PMN-PT-BT  at  5°C 


(a)  Polarization  as  function  of  (b)  Electric  field-induced  strain 


Numerical  values  of  model  parameters  for  PMN-PT-BT  at  5°C  are  given  in  Table  1.  Dependence 
of  polarization  and  induced  strain  on  electric  field  is  presented  in  Figure  7. 

The  above  model  is  valid  for  small  strains,  so  consequently  we  consider  a  problem  where  induced 
strains  are  small,  in  which  case  stress,  strain,  electric  displacement  and  electric  field  coincide  in 
initial  and  deformed  configurations. 


7.1  Actuator  example 

This  example  demonstrates  the  ability  of  mathematical  homogenization  to  model  the  dependence 
of  coarse-scale  behavior  on  fine-scale  details. 


Consider  a  beam  with  top  half  made  of  electroactive  material  and  lower  half  made  of  material  with 
no  electromechanical  coupling  (Figure  8).  The  electroactive  material  is  a  periodic  composite  with 
unit  cell  consisting  of  matrix  material  and  horizontal  electroactive  fiber.  Both  materials  have  the 
same  mechanical  properties  (linear  isotropic  material  with  Young  modulus  E  and  Poisson  ratio  v 
),  with  fiber  additionally  having  electromechanical  coupling  -  electrostriction  with  parameters 
given  in  Table  1.  Unit  cell  dimensions  are  100x100. 


The  left  side  of  the  beam  is  mechanically  fixed  and  grounded.  A  potential  36GV  is  applied  to  the 
right  hand  side  of  the  beam. 


(a)  undeformed  beam  (b)  unit  cell 


Figure  8:  a)  geometry  of  the  beam  and  boundary  conditions:  the  upper  half  is  made  of  heterogeneous 
material,  containing  electroactive  phase  and  phase  without  electromechanical  coupling  and  consisting  of 
periodic  repetition  of  unit  cells,  b)  a  unit  cell  consisting  of  matrix  and  electrostrictive  horizontal  fiber,  c) 
deformation  of  the  beam  when  voltage  is  applied  to  the  right  side 


When  voltage  is  applied,  the  electrostrictive  material  in  the  upper  half  of  the  beam  extends  due  to 
electrostriction  effect  and  the  beam  will  bend  as  shown  on  figure  8.c.  This  permits  using  the  beam 
as  an  actuator. 


We  considered  three  cases  with  different  fiber  thicknesses  (diameter  in  3D)  a:  a=100  (i.e.  the 
whole  unit  cell  consists  of  electrostrictive  material),  a= 60  and  a=40.  The  deflections  of  the  point  A 
for  different  fiber  thicknesses  a  found  by  DNS-simulation  are  shown  in  Figure  9. a.  It  can  be  seen 
that  different  fiber  sizes  result  in  different  actuation  capabilities  for  the  same  loading,  i.e.,  coarse- 
scale  response  is  sensitive  to  the  fine-scale  details. 

(a)  (b) 


(d) 


Figure  9:  (a)  displacements  of  point  A,  DNS  simulations  of  actuation  with  different  fractions  of 
electrostrictive  phase  and  Comparison  of  DNS  (red  line)  and  MH  (blue  line  for  96  element  and  green  line 
for  192  elements)  solutions,  for  different  thicknesses  of  electroactive  fiber  b)  o=100,  c)o=40  d)  a= 20. 


The  purpose  of  homogenization  is  to  capture  this  dependence.  Figure  9  shows  also  the  comparison 
of  DNS  simulation  with  MH  for  different  values  of  a.  The  MH  simulations  were  performed  for  two 
different  meshes  with  96  elements  and  with  192  elements  to  study  the  convergence  of  MH  solution 
to  DNS  solution.  For  a=20,  the  error  between  DNS  and  MH  was  4.8%  for  96  element  mesh  and 
1.3%  for  192  element  mesh.  For  a=40,  the  error  was  4.9%  and  1.4%,  respectively.  For  a=100, 
error  was  5%.  These  results  suggest  that  the  mathematical  homogenization  is  capable  of 
reproducing  the  reference  solution  and  captures  the  dependence  of  the  coarse-scale  response  on 
fine-scale  details. 

7.2  Beam  bending 

Now  consider  the  whole  beam  made  of  an  electrostrictive  material  subjected  to  a  linearly  varying 
traction  F  =  0.25-ldE'  +  5  ■  TPa,  as  shown  in  Figure  10.  The  left  side  of  the  beam  is  mechanically 


fixed  and  grounded.  Electric  potential  cp  is  applied  to  the  right  hand  side  of  the  beam.  Again,  the 
unit  cell  consists  of  matrix  and  fiber  as  depicted  in  Figure  8.b.  with  a=40.  The  matrix  is  assumed  to 
be  hyperelastic  material  with  stress  depending  exponentially  on  the  first  invariant  of  the  strain 

A  “  £U  +£22  • 


Figure  10.  Beam  with  mechanical  load  applied  to  the  left  end. 

The  application  of  electric  field  changes  the  mechanical  properties  of  the  beam  which  results  in 
different  deflections  under  the  same  loading.  We  compare  the  results  of  DNS  and  MH  to  show  the 
ability  of  MH  to  capture  the  change  of  mechanical  properties  due  to  biasing  fields. 

Figure  1 1  shows  the  deflection  of  point  A  for  the  two  cases.  In  the  first  case  the  applied  potential 
was  r/7=36GV.  The  result  of  DNS  simulation  is  depicted  by  red  line  and  MH  (192  elements) 
simulation  is  depicted  with  a  green  line.  The  difference  between  DNS  and  MH  at  time  t=0.5  was 
5.6%  for  96  elements  and  1.6%  for  192  elements.  In  the  second  case  the  right  end  of  the  beam  was 
grounded,  (p=  0.  The  result  of  DNS  simulation  is  depicted  by  blue  line  and  MH  (192  elements) 
simulation  is  depicted  by  magenta  line.  The  difference  between  DNS  and  MH  simulation  at  time 
t=0.5  was  5.2%  for  96  elements  and  1.5%  for  192  elements. 

It  can  be  seen  that  the  mathematical  homogenization  is  capable  of  capturing  the  change  of 
mechanical  properties  due  to  biasing  electric  fields  and  reproduces  the  reference  solution  with 
good  approximation. 


Figure  11.  Comparison  the  DNS  and  MH  simulations  for  beam  bending. 


8.  Conclusions  and  future  work 

The  examples  considered  in  this  manuscript  show  that  nonlinear  mathematical  homogenization 
captures  well  the  coarse-scale  behavior  of  heterogeneous  electroactive  composite  and  its 
dependence  on  the  fine-scale  details. 

At  the  same  time  the  method  shares  the  shortcomings  common  to  the  first  order  homogenization 
methods;  it  is  insensitive  to  the  absolute  size  of  the  UC  because  of  assumption  of  UC 
infinitesimality  [81].  This  lack  of  accuracy  increases  with  the  unit  cell  size  and  the  magnitude  of 
strains  and  electric  fields  inhomogeneities. 

In  future  we  plan  to  study  mathematical  homogenization  for  various  dynamic  problems.  In 
particular  of  interest  are  studies  of  wave  propagation  in  the  media  with  periodic  resonant 
structures.  It  would  be  interesting  to  explore  if  mathematical  homogenization  can  capture  the 
negative  effective  refractive  indexes  as  it  was  found  in  metamaterials.  Can  this  framework  be  used 
to  control  bandgaps  by  biasing  fields  and  will  it  allow  extending  effective  bandgaps  for  use  in 
various  devices  (for  subwavelength  imaging,  wave  attenuation  etc).  The  other  interesting 
phenomenon,  which  possibly  may  be  studied  using  mathematical  homogenization  is  the  controlled 
response  of  smart  structures  to  various  impacts.  Depending  on  the  impact  a  control  strategy  may  be 
developed  to  obtain  the  desired  response  from  the  smart  structure.  This  phenomenon  might  be 
utilized  in  development  of  smart  armor  and  other  structures.  Fine-scale  structure  optimization 
aimed  at  optimizing  coarse-scale  properties,  in  particular  for  electroactive  polymers  containing 
composites,  is  another  useful  application  of  the  method.  The  asymptotic  expansion  would  allow  to 


isolate  the  term  (p  |x,  Y,/j ,  responsible  for  electrostrictive  and  electrostatic  coupling.  Making  this 
term  larger  will  permit  increasing  the  actuation. 
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